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ABSTRACT 

We present an analysis of "bistability" in gas-phase chemical models of dark 
interstellar clouds. We identify the chemical mechanisms that allow high- and 
low-ionization solutions to the chemical rate-equations to coexist. We derive 
simple analytic scaling relations for the gas densities and ionization rates for 
which the chemistry becomes bistable. We explain why bistability is sensitive 
to the Hjj" dissociative recombination rate coefficient, and why it is damped by 
gas-grain neutralization. 

Subject headings: ISM:molecules - molecular processes 

1. Introduction 

In this paper we reexamine and analyze the phenomenon of "bistability" that occurs 
in chemical models of dark interstellar clouds. Bistability was first discussed by Le Bourlot 
et al. (1993), who found that multiple solutions to the chemical rate-equations sometimes 
appear in numerical computations of the gas-phase abundances of atomic and molecular 
species for steady-state dark cloud conditions. When bistability occurs, three sets, rather 
than a single unique set, of chemical abundances are predicted for identical input cloud 
parameters such as the total gas density, the cosmic-ray (or X-ray) ionization rate, and the 
gas-phase supply of heavy elements. Under such circumstances, the solution that is converged 
to in a numerical simulation, and perhaps, the actual steady-state of a real physical system, 
depend on the initial conditions and history of the gas. 

The subject of bistability begins with Pineau des Forets et al. (1992) who studied how 
the computed chemical states of dark clouds depend on the gas density and ionization rate. 
They found, as would be expected, that for a given ionization rate the fractional ionization 
of the gas decreases roughly as the square-root of the gas density. However, they also 
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found that the ionization fraction and associated chemical composition changes abruptly, or 
quasi-discontinuously, at a specific critical transition density. Below the critical density the 
clouds are in what Pineau des Forets et al. referred to as a "high- ionization" phase (HIP), 
characterized by high abundances of atomic ions. For densities greater than the critical 
density the clouds are in a "low- ionization" phase (LIP), with high abundances of molecular 
ions. Pineau des Forets et al. also identified a chemical instability that, they argued, is 
responsible for the abrupt transition that they encountered in their calculations. Bistability 
was subsequently discovered by Le Bourlot et al. (1993) who found that instead of a phase 
change at a specific critical density, the LIP and HIP "branches" can sometimes overlap for 
a narrow range of densities. A third (unstable) branch also appears, intermediate between 
the LIP and HIP solutions. 

Further investigations have explored the effects of the assumed gas phase elemental 
abundances, ionization rates, gas temperatures, chemical reaction rate coefficients, and gas- 
grain interactions, on the existence and character of the multiple solutions (Le Bourlot et 
al. 1995ab; Shalabiea & Greenberg 1995; Lee et al. 1998; Pineau des Forets & Roueff 2000; 
Viti et al. 2001; Charnley & Markwick 2003). These studies have shown that bistability is 
indeed related to the transition from the HIP to LIP (or vice versa), and that it is sensitive 
to quantities such as the gas-phase heavy-element abundances (Le Bourlot et al. 1995; Viti 
et al. 2001), and the adopted chemical network (Lee et al. 1998), including specifically the 
dissociative recombination rate coefficient of H^ (Pineau des Forets & Roueff 2000). In this 
paper we explain why the HIP and LIP can overlap, and we identify the chemical mechanisms 
that allow the possibility of multiple solutions. 

We will show that bistability can be simply understood and analyzed using just a hand- 
ful of reactions. In a nutshell, when bistability occurs atomic ions (usually S + ions) are 
the dominant positive charge carriers for both HIP and LIP conditions. Two modes of 
neutralization are then available to the gas. In the HIP, neutralization proceeds by slow 
radiative recombination. In the LIP, neutralization occurs by rapid dissociative recombina- 
tion, via the formation of intermediate molecular ions in reactions with O2 molecules (e.g., 
S + + O2 — > SO + + O). Bistability occurs when a low O2 abundance in the HIP can be main- 
tained up to maximal gas densities that are greater than the minimal densities above which 
a large 2 density is maintained in the LIP. We will show by means of a simple analysis why 
this is possible. The explanation involves a chemical instability, of a type first introduced 
by Pineau des Forets et al. (1992), but not given sufficient attention in subsequent studies. 

Bistability may be of primarily theoretical interest since, as we will discuss, it is a 
purely gas-phase chemical effect. If realistic gas-dust neutralization processes are included, 
the multiple solutions are damped or eliminated entirely (cf. Shalabiea & Greenberg 1995). 
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We will illustrate this by considering the effect of grain-assisted ion-electron recombination 
(Weingartner & Draine 2001) on the chemistry. Nevertheless, because bistability has been 
the subject of considerable discussion in the literature, we feel that a fresh analysis of how, 
why, and when it occurs is warranted. 

The plan of our paper is as follows. In §2 we present results of an illustrative numerical 
computation in which bistability occurs. In §3 we discuss the key chemical instability that 
drives the steady-state solutions to either the LIP or HIP. In §4 we derive simple analytic 
formulae that show how the range of gas densities and ionization rates for which the gas 
becomes bistable depends on the elemental abundances and several important chemical rate 
coefficients. In §5 we analyze the sensitivity of bistability to the dissociative recombi- 
nation rate coefficient. In §6 we explain why gas-grain neutralization processes damp and 
eliminate bistability. We summarize in §7. 



2. Computations 

We first illustrate the phenomenon of bistability in a detailed numerical computation. 
For this purpose we calculate the steady-state atomic and molecular abundances for dark 
cloud conditions, as functions of the gas density. We include only gas-phase formation 
and destruction processes, except for molecular hydrogen (H 2 ), which we assume is formed 
on grain surfaces. The chemistry is driven by cosmic-ray impact ionization of hydrogen and 
helium, followed by standard sequences of ion-molecule, charge-exchange, and neutral-neutral 
reactions. The chemistry is also mediated by radiative and dissociative recombination, and 
cosmic-ray induced photoionization and photodissociation. The effects of externally incident 
FUV photons are excluded. 

We use the Boger & Sternberg (2005) network for a system consisting of 70 atomic and 
molecular carbon, nitrogen, oxygen, sulfur, and silicon bearing species (see also Sternberg 
& Dalgarno 1995). We neglect all other heavy elements. Our set of 868 reactions and rate 
coefficients is based mainly on the UMIST99 compilation (Le Teuff et al. 2000) with some 
modifications and alterations 1 . As in Boger & Sternberg (2005) we adopt ( Oph gas-phase 
abundances, with C/H=1.32 x 10~ 4 (Cardelli et al. 1993), N/H=7.50 x 10~ 5 (Meyer et 
al. 1997), 0/H=2.84 x 10~ 4 (Meyer et al. 1998), Si/H=1.78 x lO" 6 (Cardelli et al. 1994), 
and S/H=8.30 x 10~ 6 (Lepp et al. 1988). We set He/H=0.1. We assume a constant gas 
temperature of T = 50 K for the entire range of gas densities that we consider. 



Our input list of rate coefficients is available at ftp://wise3.tau.ac.il/pub/amiel/pdr. 
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In our calculations we solve the equations of chemical equilibrium 

Y ki i l ( T ) n i n « + Y = Ui {Y k i aUl + Y & } ( 1 ) 

jl j ji j 

for the densities n« of all species % in our set. In these equations, kiji(T) are the temperature- 
dependent rate coefficients (cm 3 s -1 ) for chemical reactions between species j and / that lead 
to the production of i, and ^ are the cosmic-ray destruction rates (s _1 ) of species j with 
products %. The parameters ^ include impact ionizations of hydrogen and helium by the 
cosmic-ray protons and secondary electrons, as well as photo-destruction of heavy atomic 
and molecular species by induced UV photons. The ^ are all proportional to ( = lCT 17 C-i7 
s -1 , the cosmic-ray ionization rate of molecular hydrogen (Sternberg et al. 1987; Gredel et 
al. 1989) 

The densities rii must also satisfy the supplementary equations of element and charge 
conservation. These are, 

Y a im n i = X m n tot (2) 

i 

where ai m is the number of atoms of element m contained in species i, X m is the total 
abundance of element m relative to the total density n tot of hydrogen nuclei, and 

Y ^ = ( 3 ) 

i 

where is the charge of species i. In our computations, the negative charge is carried entirely 
by free electrons. 

H2 molecules are formed primarily on dust grains and are destroyed by cosmic-ray 
ionization. A separate equation may then be written for the H/H 2 balance, 

(>h 2 = Rn tot n R (4) 

where riu and n H2 are the atomic and molecular hydrogen densities respectively, and R = 
3 x 10~ 18 T 1 / 2 cm 3 s" 1 is the H 2 grain surface formation rate coefficient (Biham et al. 2002; 
Vidali et al. 2004). 

For fixed gas temperature and gas-phase elemental abundances, the density fractions 
Xi = rii/n tot as determined by equations (l)-(4) depend on the single parameter n tot /( (e.g., 
Lepp & Dalgarno 1996; Lee et al. 1998), as is readily seen by dividing each term in these 
equations by n^ ot . We consider a range of parameters for which almost all of the hydrogen is 
in the form of H 2 , so that n to t = 2nn 2 . We therefore present our results as functions of the 
ratio nu 2 /C-n (cm -3 ). 
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We solve equations (l)-(4) numerically via Newton-Raphson (NR) iteration. In our 
procedure we vary the density along the HIP and LIP branches and use the solutions from 
the previous density step as the "initial guess" for the NR iteration in the next density step. 
In the bistability region we carry out repeated NR iterations with random initial guesses, so 
that all three branches are found. 

In Figure 1 we plot the densities nj/C-17, and density fractions rii/n-u 2 , as functions of 
n H 2 /C-i7 ranging from 10 2 to 10 4 cm~ 3 , for several species of interest. Panels (a) and (b) 
display our results for electrons, C + , S + , H^", and a sum over all molecular ions. Panels 
(c) and (d) display our results for the neutral species C, O, CO and O2. Several important 
features are apparent in Figure 1, some of which have been noted in previous discussions. 

For our assumed elemental abundances and gas temperature, a high-ionization phase 
(HIP) is present for densities less than (riH 2 /C-i7)max = 675 cm -3 , and a low-ionization 
phase (LIP) is present for densities greater than (nH 2 /C-i7)mm = 370 cm~ 3 . The HIP and 
LIP overlap, and the gas is bistable, between these minimal and maximal densities. The 
HIP and LIP solutions are connected by a third branch that we refer to here as the MIX 
solution. 

In the HIP, the positive charge is carried entirely by atomic ions. S + is the dominant ion 
down to nu 2 /C-n — 60 cm -3 . At lower densities C + takes over, and for nu 2 /C-n < 11 cm -3 , 
H + dominates. In the LIP, the molecular ions HCO + and H 3 + are the primary positive 
charge carriers for densities greater than n H2 /C-i7 > 6 x 10 3 cm -3 . However, S + dominates 
at lower densities in the LIP. Most importantly, in our calculation S + is the primary positive 
charge carrier in both the HIP and LIP for the densities at which bistability occurs. 

The fractional ionization increases from a minimum value x e = 6 x 10~ 6 as the density 
decreases from (riH 2 /C-i7)max i n the HIP. In the LIP, the fractional ionization decreases from 
a maximum value of x e = 2 x 10~ 6 as the density increases from (riH 2 /C-i7)min- The fractional 
ionization drops by a factor of ~ 10 in moving from the HIP to LIP in the bistability region. 

The free atomic carbon density and the C/CO abundance ratio is large in the HIP. The 
ratio C/CO 0.3 in the bistability region, and exceeds unity for n H2 /C-i7 ^ 55 cm -3 . In 
the LIP the C density and the C/CO ratio are small. In our computation the C/CO ratio 
drops by a factor ~ 50 in moving from the HIP to the LIP. The drop in C/CO is associated 
with an increase in the O2 density and the O2/CO ratio. It is very significant that the C 
and 2 densities vary in opposite directions in the transition between the HIP and LIP. 

On the MIX branch the fractional ionization, and the C/CO and O2/CO ratios, are 
intermediate between the HIP and LIP. S + is also the dominant positive ion on the MIX 
branch. However, the fractional ionization increases with gas density along the MIX. 
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In Table 1 we list the densities of selected atomic and molecular species for the three 
HIP, LIP, and MIX solutions, for n H2 = 4.3 x 10 2 cm" 3 and (=1.0x 1(T 17 s _1 . We also 
list the HIP/LIP density ratios for these values of % 2 and (. Atomic ions, e.g., S + , Si + , 
and C + , are abundant in the HIP solution. In the LIP solution, molecular ions such as Hj~, 
HCO + , and H30 + are abundant, leading to more OH, H2O, and O2 than in the HIP. The 
oxygen bearing molecules, NO, SO, and SiO, are enhanced by the large O2 density that is 
maintained in the LIP. Carbon bearing molecules, such as CN and HCN, are enhanced in the 
HIP because of the large atomic carbon density in this phase (see also Boger & Sternberg 
2005). CS is produced by reactions of SO (enhanced in the LIP) with C (enhanced in the 
HIP) so that the CS density changes by only a factor of ~ 2 in moving from the LIP and 
HIP. 



3. Chemistry and Instability 

Figure 1 shows that as the gas density n H2 is increased along the sequence of HIP 
solutions an abrupt jump to the LIP branch occurs when % 2 rises above (wH 2 /C-i7)max- 
Similarly, as nn 2 is decreased along the LIP, a jump to the HIP occurs when nn 2 falls 
below (nH 2 /C-i7)min- As we now discuss, these jumps are signatures of a simple chemical 
instability involving electrons, S + and ions, C atoms, and O2 molecules. The instability 
can be understood by considering the gas-phase formation and destruction sequences for 
these species, for gas densities where the bistable solutions occur. We first discuss the HIP 
and LIP, and we then consider the MIX. 

We begin with H3". In both the HIP and LIP, this ion is formed by cosmic-ray ionization 
of H 2 followed by proton transfer, 

H 2 + cr^H+ + e (Rl) 

H+ + H 2 ^H+ + H . (R2) 
In the HIP, the electron fraction is large, and dissociative recombination 

H+ + e -> H 2 + H (R3) 

H+ + e^H + H + H (R4) 

dominates the removal of the H3 . However, in the LIP the electron fraction is small, and 
the is removed by 

H+ + CO -> HCO+ + H 2 . (R5) 
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For a gas temperature of 50 K, dissociative recombination dominates the destruction of 
when n(e)/ra(CO) > 10" 2 . 

In both the HIP and LIP, the production of 2 is initiated by proton transfer 

O + H+ -> OH+ + H 2 , (#6) 

which is followed by abstraction 

OH+ + H 2 -»■ H 2 0+ + H (R7) 
H 2 0+ + H 2 -> H 3 0+ + H , (R8) 

dissociative recombination 

H 3 0+ + e ->• OH + H 2 (#9) 
H 3 0+ + e^H 2 + H , (RIO) 

and the neutral-neutral reaction 

OH + 0^0 2 + H . (mi) 

However, the 2 is destroyed in different ways in the two phases. In the HIP, the density of 
free atomic carbon is large, and the 2 is removed by 

2 + C -> CO + O . (#12) 

In the LIP, the C density is small, and the 2 is removed much less efficiently in reactions 
with less abundant atoms and atomic ions. An upper limit on the 2 density in the LIP is 
set by cosmic-ray induced photodissociation 

2 + crp -> O + O . (#13) 

Reaction (R13) is inefficient, and the 2 density becomes large in the LIP. 

In the HIP, C + is produced by cosmic-ray driven helium impact ionization of CO, 

He + cr -> He + + e (RU) 

He + + CO -> C + + O + He . (#15) 
This is then followed by rapid charge transfer 

S + C + -> S + + C (#16) 
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which leads to the production of S + and C. The S + ions are removed by radiative recombi- 
nation, 

S+ + c -> S + v . (#17) 
The C atoms are removed by cosmic-ray induced photoionization 

C + crp -> C+ + e . (#18) 

Because cosmic-ray photoionization is inefficient, a large C density is maintained in the HIP. 

In the LIP, C + is also produced by helium impact ionization ([R14] and [R15]). However, 
because the 2 density is large in the LIP, the C + ions are rapidly neutralized by the sequence 

C + + 2 -> CO + + O (#19) 

CO+ + e -> C + O , (#20) 

and the C + density becomes small. In the LIP, charge transfer between C + and S ([R16]) 
contributes negligibly to the removal of C + , or to the formation of S + . The main source of 
S + in the LIP is cosmic-ray induced photoionization 

S + crp^S + + e . (#21) 

Crucially, 2 also dominates the neutralization of S + in the LIP, via 

S+ + 2 -> SO+ + O (#22) 

SO+ + e -> S + O . (#23) 

For a gas temperature of 50 K, neutralization of S + by 2 dominates over radiative recom- 
bination when n(0 2 )/n(e) > 0.8. Finally, the C atoms formed by (R20) are also removed by 
reactions with 2 ([R12]), and the C density becomes small in the LIP. 

For the gas densities of interest, most of the gas-phase sulfur is present as S and S + in 
both the HIP and LIP, with only a small fraction bound in molecules. S + is the primary 
positive charge carrier in both the HIP and LIP. However, in the HIP S + is removed by 
radiative recombination, whereas in the LIP by reactions with 2 . In the HIP, H3 is removed 
by dissociative recombination, C is removed by cosmic-ray induced photoionization, and 2 
is removed by C. In the LIP, H3 is removed by CO, C is removed by 2 , and 2 is limited 
by cosmic-ray induced photodissociation. The abundance ratio C/0 2 is large in the HIP, 
and is small in the LIP. 

Now consider a gas where S + is again the dominant charge carrier, but which differs from 
the HIP or LIP in that dissociative recombination contributes significantly to the removal 
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of H3 (as in the HIP), while on the other hand S + is neutralized by 2 (as in the LIP). 
Assume also that 2 is produced by the sequence (R6)-(R11), and that the 2 destruction 
rate increases if the 2 density decreases. This is possible when reaction (R12) between 2 
and C is the dominant destruction mechanism for both species, so that n(C) increases if 
n(0 2 ) decreases. For these conditions a chemical instability can occur, as we demonstrate 
with the following perturbation analysis. 

Consider a possible solution to the steady-state rate equations. With the above assump- 
tions this requires n(e) = n(S + ). Now consider a perturbative increase in the ionization 
fraction x e = n(e)/nn 2 , carried out at fixed gas density nn 2 . This will reduce by a factor 
that depends on the relative efficiencies of dissociative recombination versus reaction with 
CO in removing the H3 . The 2 is removed by C (via [R12]), so for a given C density the 
reduction in Hjj" will lead to an equal reduction in 2 (since 2 is formed by [R6]-[R11]) 
2 . However, the decrease in 2 will lead to an increase in C, since the C is removed by 
2 (again via [R12]). This reduces 2 further still. The S + density then increases by a 
factor equal to the total reduction in 2 (since S + is removed by 2 ). Finally, to close the 
loop, n(e) must then be set equal to the new value of n(S + ), since by assumption S + carries 
the positive charge. However, if the increase in S + is larger than the initial increase in the 
electron density the cycle "runs away" in a positive feedback loop. 

This H 3 f -0 2 -S + cycle is illustrated schematically in Figure 2. It is clear that the cycle can 
also lead to a run-away decrease in the electron density, in response to a perturbative decrease 
in n(e). For increasing electron densities, 2 decreases and C increases. For decreasing 
electron densities, 2 increases and C decreases. 

The instability shuts off for sufficiently large or small fractional ionizations. For large 
x e , radiative recombination, rather than removal by 2 , becomes the dominant S + neutral- 
ization process. For small x e , reaction with CO replaces dissociative recombination as the 
Hjj" removal mechanism. Either way, the H 3 f -0 2 -S + instability cycle shuts off, and the gas 
becomes chemically stable. 

In the HIP, dissociative recombination dominates the destruction of H3", but S + is 
removed by radiative recombination. In the LIP, S + is removed by 2 , but H^ is removed by 
CO. Crucially therefore, the H3-0 2 -S + cycle is inoperative on the HIP and LIP branches, 
and these solutions are chemically stable. 

We now identify the maximum possible density for the HIP, (nH 2 /C-i7)ma X ) a s the point 
at which 2 , as formed and destroyed in the HIP, becomes sufficiently abundant to begin 



2 In this analysis we assume that most of the oxygen not locked in CO remains atomic. 
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removing the S + . When this occurs, the Hjj"-0 2 -S + cycle turns on (since H3 is still being 
removed by electrons) and a stable HIP is no longer possible. As the density is increased 
above (riH 2 /C-i7)max the instability induces a sharp drop in 71(e). As the electron density 
drops in the run-away loop, O2 increases and C decreases. The system is finally stabilized 
when the electron density becomes sufficiently small for reactions with CO to begin domi- 
nating the removal of H^ . The H^-C>2-S + cycle then turns off, and a stable solution on the 
LIP is reached. 

Similarly, we identify the minimum possible density for the LIP, (nH 2 /C-i7)mm, as the 
point at which the fractional ionization, as determined for the LIP, becomes sufficiently large 
for dissociative recombination to begin dominating the removal of H^. When this occurs, 
the H3 -0 2 -S + cycle turns on, (since S + is still removed by 2 ) and a stable LIP is no longer 
possible. As the density is decreased below (ra^/C-r^min the instability induces a sharp 
rise in 71(e). As the electron density increases, 2 decreases and C increases. The electron 
density rises until it becomes sufficiently large (and the 2 density sufficiently small) for 
radiative recombination to dominate the neutralization of S + . The instability then shuts off, 
and a stable solution on the HIP is reached. 

The MIX branch can now also be understood. It represents the locus of equilibrium 
solutions, with varying tih 2 , for which the H 3 f -0 2 -S + cycle is operative, but for which any 
change in electron density is precisely matched by a corresponding change in the S + density. 
For these solutions, dissociative recombination and reactions with CO both contribute to 
the removal of H^ , and S + is removed by 2 . Furthermore, reaction (R12) is the dominant 
destruction mechanism, for both 2 and C. As the electron fraction increases along the MIX, 
the 2 formation rate decreases, but by a factor that is smaller than the increase in electron 
density. The density of C atoms (which remove the 2 ) increases such that the total increase 
in 2 , and the resulting increase in S + , exactly matches the increase in electron density. 

Figure 3 summarizes the key characteristics of the HIP, LIP, and MIX branches, as 
described above. The triangular symbols indicate the three "legs" of the H3 -0 2 -S + cycle 
illustrated in Figure 2. In the HIP, removal of S + by 2 is inoperative so this leg is labeled 
as "OFF" . In the LIP, dissociative recombination of H3 is inoperative so this leg is "OFF" 
for this phase. For the MIX, which shares properties of both the HIP and LIP, all three legs 
of the H 3 f -0 2 -S + cycle operate, and are "ON". In Figure 3 we also indicate the dominant 
removal reactions for S + , H3 , 2 , and C, for each of the three branches. 

The HIP and LIP as originally defined by Pineau des Forets et al. (1992), were dis- 
tinguished by whether H^ is removed by dissociative recombination (HIP), or by reactions 
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with CO (LIP) 3 . Here we have identified a second, and crucial, distinction. In the HIP, 
atomic ions, are removed by (slow) radiative recombination. In the LIP the atomic ions are 
effectively removed by rapid dissociative recombination, via the formation of intermediate 
molecular ions, as in the sequence (R22) and (R23). 

In our discussion it is S + that plays the crucial role of positive charge carrier in both the 
HIP and LIP. The reason why simultaneous HIP, LIP, and MIX solutions are possible now 
emerges. Bistability occurs when radiative recombination is able to dominate the removal 
of S + in the HIP up to maximal gas densities that are greater than the minimal densities 
down to which reactions with CO dominate the removal of Hjj" in the LIP. Alternatively, 
bistability occurs when it is possible for a low HIP O2 abundance to be maintained up to 
maximal gas densities that are greater than the minimal gas densities down to which a large 
LIP 2 density can be maintained. This is our key result. When this occurs, there exists a 
range of densities for which the S + ions are either removed slowly via radiative recombination 
(HIP, low O2) or rapidly via dissociative recombination (LIP, high O2). The simultaneous 
solutions are separated by the Hj~-02-S + instability cycle, which drives the O2 density to 
either small (HIP) or large (LIP) values. 

The critical densities (riH 2 /C-i7)mm and (riH 2 /C-i7)max, as defined above, depend on a 
variety of factors including the elemental abundances, and several specific chemical rate 
coefficients. In §4 we derive simple analytic expressions for these critical densities. We 
determine, analytically, when it is possible for (riH 2 /C-i7)max to exceed (?iH 2 /C-i7)mm, which is 
our condition for bistability, and we determine the range of gas densities for which bistability 
occurs. 



4. Analysis 

4.1. HIP and n max 

For stable HIP conditions the electron density must be sufficiently large for dissociative 
recombination to dominate the removal of Hj~, and the 2 density must be sufficiently small 
for radiative recombination to dominate the removal of S + . We now show analytically that 
consistent solutions that satisfy these HIP conditions are possible provided the density is 
less than a critical maximum gas density. 



3 These authors also recognized that the transition from the LIP to HIP is driven by a chemical instability 
involving 02- However, they did not consider how the instability is quenched in each phase. 
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We first derive an expression for the 2 density in the low density limit of the HIP. We 
show that in this limit the O2 density is small, and radiative recombination dominates the 
removal of S + . We then show that our low density expression for O2 remains consistent only 
up to a maximal gas density. We identify this as the maximum possible density for HIP 
solutions. We then verify that dissociative recombination dominates the removal of Hj, a 
necessary condition for the HIP, for all densities below this critical density. 

In the HIP, O2 is produced by the sequence (R6)-(R11), and is removed by reaction 
(R12) with atomic carbon. This gives, 

_ n(H+)n(Q)MoH /,x 

where k$ = 8.0 x 10~ 10 and ki 2 = 3.3 x 10~ n cm 3 s _1 are the rate coefficients (at 50 K) of 
(R6) and (R12), and g?oh = 0-74 is the fraction (Jensen et al. 2000) of H 3 + recombinations 
that produce OH. Atomic carbon is produced by helium impact ionization of CO followed 
by immediate charge transfer neutralization ([R15] and [R16]), and is removed by cosmic-ray 
induced photoionization ([R18]). Because (R15) dominates the removal of He + , the rate of 
C formation equals the cosmic-ray ionization rate of helium. Therefore, 

n>H 2 pis 

independent of nn 2 /C- 111 this expression p± 8 = 1.02 x 10 3 is the efficiency of (R18), and 
X He = 0.1 is the helium abundance. It is clear that a large C density is indeed expected 
in the HIP, since n(C)/nn 2 = 10~ 4 is comparable to the typical total gas-phase carbon 
abundance (e.g. Xq = 1.32 x 10~ 4 in our numerical model). The large C density keeps the 
O2 abundance small in the HIP. 

The 2 density is proportional to Hj" , which in turn, depends on the electron density. 
For electrons we set n(e) = n(S + ). The S + is produced the same way as atomic carbon, 
i.e., by helium impact ionization of CO followed by charge transfer ([R15] and [R16]). We 
assume that the S + is removed by radiative recombination. This gives, 

, e ^=(i«) 1/2 = 2.9xlO-(^I) 1/2 . (7) 
n H2 \k 17 n H2 y Vn H2 / 

where the radiative recombination rate coefficient kn = 1.2 x 10~ n cm 3 s -1 . In this limit 
the fractional ionization varies as (^h 2 /C)~ 1//2 - F° r n H 2 /C-i7 = 100 cm -3 , equation (7) gives 
x e = 2.9 x 10 -5 in good agreement with 2.3 x 10~ 5 computed numerically (see Figure 1). 

H3 is produced by cosmic-ray ionization and is removed by dissociative recombination, 
so that 

= -L (-^) 1/2 .2.0X10-^) 1/2 ,8) 
nu 2 h A n{e) k 3A \X He nu 2 y v n H2 / 
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where fc 3j 4 = 1.7 x 1CT 7 cm 3 s _1 is the total rate coefficient for dissociative recombination. 
The Hg abundance also varies as (n^/C)' 1 ^ 2 - F° r n n 2 /C-i7 — 100 cm -3 , equation (8) gives 
n(Fhj~)/nH 2 = 2.0 x 10~ 8 , as also computed numerically. 

Inserting equations (6), (7), and (8) into equation (5) gives 



n H2 fcl2 ^ Vn H 2 / V^H : 

where Xq = Xq — Xq = 1.5 x 10~ 4 , and where we have assumed that most of the oxygen not 
bound in CO remains atomic. The O2 abundance varies as (^h 2 /C)~ 1 ^ 2 - F° r n H 2 /C-i7 = 100 
cm~ 3 , equation (9) gives 71(02) /wh 2 — 11 x 10~ 6 , close to 3.0 x 10~ 7 computed numerically. 

The O2 density as given by equation (9), is sufficiently small for removal of S + by O2 
to be negligible compared to radiative recombination. This can be seen by considering the 
ratio of the S + removal rates 



k 22 n(0 2 ) 2Xq k 6 k 22 



k 17 n(e) X^ e k 3A k 12 



(Wis = 0.05 (10) 



where k 22 = 1.5 x 10~ n cm 3 s _1 is the rate coefficient for reaction (R22). Because this ratio 
is independent of nu 2 /(, the removal of S + by O2 is negligible at all densities so long as 
equation (9) for the O2 abundance is an appropriate approximation. 

Equation (9) is valid so long as cosmic-ray photons dominate the removal of the carbon 
atoms, so that the C density is given by equation (6). However, cosmic-ray induced destruc- 
tion occurs at a rate independent of nn 2 , whereas the O2 density as given by equation (9) 

1 /2 

varies as n Ha . Therefore, for sufficiently large nn 2 reactions with O2 ([R12]) must eventually 
dominate the removal of the C atoms. When this occurs the C density becomes small, and 
the O2 density must become large (see equation [5]). We identify this as the point where O2 
becomes sufficiently abundant to begin removing the S + . However, because H3 continues to 
be removed by dissociative recombination, the Hjj"-02-S + instability cycle turns on at this 
point, and a stable HIP is no longer possible. 

Thus, the HIP can be maintained so long as the O2 density given by equation (9) is 
small enough for cosmic-ray photoionization to dominate the removal of the C atoms. This 
condition, 

hMO^ _ 2X^^ do jn^ < 1 (n) 



provides a maximal density for HIP solutions. For our standard elemental abundances and 
rate-coefficients 

= 7.63 x 10 2 cm" 3 (12) 



( n n 2 \ 
VC-17A 
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This is slightly larger than our numerically computed density of 675 cm 3 at which the 
transition to the LIP occurs (see Figure 1). 

For consistency we must verify that dissociative recombination dominates the removal 
of H3" for densities as high as (nH 2 /C)max- Dissociative recombination dominates when the 
ratio 

hMe) fXu c \^ 2 k 3A 1 / C \ 1/2 



A; 5 n(CO) \ k 17 J k 5 X c \n n J >> 1 ' (13) 

where k^ = 1.7 x 10~ 9 cm 3 s -1 is the rate coefficient of reaction (R5). For (n/C-i7) m ax = 
7.63 x 10 2 cm" 3 this ratio equals 8.0, so consistency is obtained. 



4.2. LIP and n min 

In the LIP, the electron density must be sufficiently small for CO to dominate the 
removal of ILt, and the O2 density must be sufficiently large for O2 to dominate the removal 
of S + . We now show that these conditions can be satisfied, provided the gas density is greater 
than a critical minimum density. 

To show this, we first derive an expression for the electron density assuming that S + is 
removed by abundant O2, with an upper limit for the O2 density set by cosmic-ray induced 
photodissociation. We show that CO then dominates the removal of Hjj" for densities greater 
than a minimum critical density. We then verify that the assumption that O2 dominates the 
removal of S + is consistent for densities greater than this critical density. 

For the electron density in the LIP we again assume that n(e) = n(S + ). In the LIP, 
S + is produced by cosmic-ray induced photoionization, and we assume that it is removed by 
2 . This gives, 

n(e) _ 2X s p 2 iC 



nn 2 k 2 2n{0 2 ) 

where p 2 i = 2.0 x 10 3 is the efficiency of (R21). H3 is removed by CO, so that 



(14) 



n ( H ^) _ C _ 1 L = 2 .2xl0- 5 ^ , (15) 

n Ha k 5 n(CO) 2k 5 X c n n2 n Ha 

where we have assumed that in the LIP all of the carbon is locked in CO. For nn 2 /C-i7 = 10 3 > 
equation (14) gives « H +/riH 2 = 2.2 x 10~ 8 , in good agreement with the numerically computed 
1.0 x 10~ 8 . 



For 2 we again assume the production sequence (R6)-(R11). If the 2 is removed by 
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only cosmic-ray induced photodissociation, 



n(0 2 ) n(H+)n(0)MoH h 1 J 00win _ 4 



XT — c^oh^ = 2 - 3 x 10 > ( 16 ) 

™H 2 CPl3™H 2 I C 

where p 13 = 1.7 x 10 3 is the efficiency of (R13). This is actually an upper-limit for the 
O2 abundance because additional chemical reactions contribute to the removal of the O2 in 
the LIP. In this approximation the O2 density is independent of nu 2 /(, consistent with our 
numerical computations (see Figure 1). Because the upper limit is close to the total available 
oxygen not locked in CO, we set n(0 2 )/n H2 = 2/o 2 ^o? where the 2 fraction f Q2 < 1. In 
our numerical computation f Q2 « 0.2 (see Figure 1). 

The electron fraction, x e , is then 

It = wtt- f = 1A x 10 "f it (17) 

n H2 A Q /C22 Jo 2 n H 2 /o 2 n H 2 

In this limit x e varies as (^h 2 /C) _1 - For f 02 = 0.2, equation (17) gives x c = 3.7 x 10~ 7 , in 
good agreement with 4.1 x 10~ 7 in Figure 1. 

The electron density as given by equation (17) is independent of % 2 , whereas the CO 
density decreases with hh 2 . Reactions with CO can therefore dominate the removal of 
only down to a minimal density set by the condition 

hn(CO) 2X C X' k 5 k 22 n U2 

1 TT = — v 1 Jo 2 —- > 1 • l iy ) 

This relation provides a minimal critical density 

(£5L) = ii cm- (19) 

\(,_17/min Jo 2 

for which LIP solutions are possible. For fo 2 = 0.2, this is slightly smaller than the numeri- 
cally computed value of 370 cm~ 3 at which the transition from the LIP to HIP occurs (see 
Figure 1). When the gas density falls below the critical density electrons remove FLi , but 
2 continue to remove S + . The H3-0 2 -S + instability cycle therefore turns on, and the gas 
is driven from the LIP to the HIP. 

For consistency we must now verify that reactions with O2 dominate the removal of S + 
everywhere in the LIP. This requires that 

fe 2 2 n(Q 2 ) _ 2X^2 k\ 2 2 n H2 2 ™h 2 

~hMey ~ ^k^ fo ~ ~ /oa c^ • (20) 

for all LIP densities, down to the minimal density. Inserting (n/C-i7) m in as given by equa- 
tion (19) gives k 2 2n(0 2 ) / ki 7 n(e) = 142/o 2 . We conclude that consistent LIP solutions are 
obtained even for / 0a as small as ~ 0.01 (for which the 2 density in the LIP is still much 
larger than in the HIP.) 
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5. Scaling 

The key result of the above analysis is that for our standard parameters (n/C-i7)max, 
as given by equation (12), is greater than (n/C-i7)min, as given by equation (19). The gas is 
therefore bistable between the maximal HIP density and minimum LIP density. The MIX 
branch must therefore "switch back" from the maximal HIP density to the minimal LIP 
density. The H 3 -0 2 -S + instability cycle drives the 2 density to either very low values in 
the HIP where a large ionization fraction is maintained by slow radiative recombination, or 
it drives the 2 density to very large values in the LIP where a low ionization fraction is 
maintained by rapid dissociative recombination. Our analytic estimates for the two critical 
densities are in good agreement with the numerical results displayed in Figure 1. 

Equations (11) and (18) show that the critical densities depend on several parameters, 
including the gas-phase elemental abundances of carbon, oxygen, sulfur and helium, and on 
several specific chemical rate coefficients, including /c 3i4 , the dissociative recombination rate 
coefficient of H^. Varying any of these parameters alters the range of densities for which 
bistability occurs, or eliminates it entirely if (n/£„i 7 ) max becomes smaller than (n/£„i 7 ) min . 

A dependence of bistability on ^3,4, a rate coefficient that has been the subject of 
much investigation and controversy (see e.g., McCall et al. 2004), was found in numerical 
computations by Pineau des Forets & Roueff (2000). In their calculations (see their Fig. 4) 
they set nn 2 /C-i7 — 2.5 x 10 3 cm -3 , and varied fc 3) 4 from 10~ 8 to 10~ 6 cm 3 s~ x . They found 
that the gas is in the LIP or HIP, for low or high values of fc 3) 4, and that it is bistable for 
intermediate Hjj" recombination efficiencies. Our analysis provides an explanation for this 
behavior. 

In Figure 4 we plot x e versus tih 2 /C-i7, f° r fi ye values of fc 3j4 , ranging from 0.2 to 5 times 
our standard 50 K value of 1.7 x 10~ 7 cm 3 s _1 . In these numerical computations, all other 
parameters, including all rate coefficients and gas-phase elemental abundances, are the same 
as in the computations we presented in §2. Figure 4 shows that for small fc 3) 4 bistability is 
suppressed, and that as A; 3 4 is increased the gas becomes bistable for an increasing range of 
densities. For nn 2 /C-i7 held fixed at, say, 1.0 x 10 3 cm -3 , the gas is in the LIP or HIP for 
low or high fc 3i4 , and is bistable for intermediate values. This is the same behavior found by 
Pineau des Forets & Roueff (2000). 

These results can be readily understood as follows. As the dissociative recombination 
rate coefficient of H^ is increased, less O2 is produced by the sequence (R6)-(R11) at any give 
density. An HIP can therefore be maintained to higher densities before O2 becomes signif- 
icant in removing S + . An increased H^ recombination efficiency also implies that electrons 
become competitive with CO in removing H^ at larger densities in the LIP. Quantitatively, 
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equation (11) shows that (nH 2 /C-i7)max, the density above which O2 begins removing S + , 
is proportional to k% 4 . However, (n H2 /C-i7)min, the density below which electrons begin re- 
moving the H^, varies only linearly with fc 3 4. Therefore, as £^4 is reduced, (ji/C-i7)max falls 
below (^H 2 /C-i7)min, and bistability disappears. As A; 3 4 is increased, (^H 2 /C-i7)max grows 
more rapidly, and eventually exceeds, (riH 2 /C-i7)min, and the gas becomes bistable for an 
increasing range of densities. For example, if k^4 is increased by a factor of 5 relative to our 
standard value, Figure 4 shows that (^H 2 /C-i7)max increases by about a factor of 25, and 
(^H 2 /C-i7)min increases by about a factor of 5. This is consistent with expressions (11) and 
(18). 

A similar scaling analysis may be carried out for the other parameters that control 
(^H 2 /C)min and (^H 2 /C)max- For example, if the oxygen abundance is reduced our analysis 
shows that the regime of bistability should move to higher values of nn 2 /(, as found nu- 
merically by Viti et al. (2001) (see their Figures 5 and 7). For a lower oxygen abundance 
( n H 2 /C)min is increased because the HIP can be maintained to a higher density before 2 
becomes abundant enough to remove the atomic ions. The reduced O2 density leads to a 
higher electron density in the LIP, so (riH 2 /C)min is increased as well. 

Our numerical and analytic results show how the onset of the H3 -0 2 -S + instability cycle 
depends on the 2 and electron densities as determined by the elemental abundances and by 
the Hjj~ dissociative recombination rate coefficient. However, Figure 4 reveals an additional 
interesting feature. For /c 3 4 equal to 0.2 times our standard value the gas again becomes 
marginally bistable at (tih 2 /C-i7) ~ 9 cm~ 3 . At such low densities the dominant atomic ion 
in the HIP is C + rather than S + (see Fig. 1). The instability, and the transition from HIP 
to LIP, is therefore modified by the addition of an identical H3 -0 2 -C + cycle, with reactions 
(R19) and (R20) operating to neutralize the C + . The behavior in this numerical example 
is a bit more complicated because S + is still the dominant ion in the LIP at the transition 
density. We note that at such low densities realistic clouds would likely become "diffuse" 
and optically thin to external far-ultraviolet radiation, with a resulting reduction in 2 and 
further modifications to the chemistry. We do not consider such effects here. 



6. Gas-Grain Recombination 

We have shown that bistability arises in purely gas-phase models, consistent with previ- 
ous findings. When the HIP and LIP solutions "overlap" , neutralization of atomic ions (S + 
in our models) proceeds by either slow radiative recombination in the HIP, or via the for- 
mation of intermediate molecular ions (S + + O2 — > SO + + O) followed by rapid dissociative 
recombination in the LIP. The two solutions are separated by a chemical instability, which 
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we have identified in our models as the Hg -0 2 -S + cycle. 

Our analysis suggests that bistability will be damped, or eliminated entirely, if additional 
neutralization processes are included that interrupt the H 3 f -0 2 -S + cycle. We demonstrate 
this by considering how the single process of "grain assisted recombination" of atomic ions 
affects the chemistry. In this process, atomic ions recombine via electron transfer from grains 
to the ions during collisions with the grains (Draine & Sutin 1987; Weingartner & Draine 
2001). The grain assisted recombination rate may be expressed as k g n tot , where n to t is the 
total density of hydrogen nuclei. For any ion, k g is an effective rate coefficient that depends on 
the dust-to-gas ratio, grain size distribution, grain charge, and gas temperature. This process 
dominates radiative recombination when 2k g > k r x e , where k r is the radiative recombination 
rate coefficient. For k r ~ 10 -11 cm 3 s _1 (appropriate for S + ), and a fractional ionization 
x e ~ 10~ 5 , as occurs for densities near (riH 2 /C-i7)max in our computation (see Fig. 1), grain 
assisted recombination dominates radiative recombination for k g > 5 x 10~ 17 cm 3 s~ 4 . This 
is well below the maximum values of k g ~ 10~ 14 cm 3 s^ 1 found by Weingartner & Draine, 
so this process is potentially relevant in determining the ionization balance. Furthermore, 
grain assisted recombination will interfere with the H 3 f -0"2-S + cycle and render it inoperative 
when kg > where as previously defined, fo 2 X' Q is the fraction of oxygen bound 

in O2, and k 2 2 = 1-5 x 10~ n cm 3 s~ 4 is the rate coefficient of reaction (R22). For /o 2 ~ 0.2 
and X' Q = 1.5 x 10~ 4 , the H 3 f -0 2 -S + instability cycle is quenched when k d >5x 10~ 16 cm 3 
s~ 4 . Bistability and the sharp transition from HIP to LIP should then disappear. 

We illustrate this in Figure 5, where we plot x e versus n^/C-n, for k g ranging from zero 
(as in Figures 1 and 4) to 1.0 x 10 -14 cm 3 s~ 4 , for all atomic ions. All other parameters are 
as in the calculation in §2. As k g is increased the diminished electron density leads to more 
H^ and 2 in the HIP, so that (nu 2 /C-n)max is reduced. Similarly, the onset of dissociative 
recombination in the LIP occurs at lower densities, and (^H 2 /C-i7)min is reduced. Thus, at 
first the bistability region shifts to lower densities as k g is increased. However, as estimated 
above, when k g exceeds ~ 3 x 10~ 16 cm 3 s" 1 bistability is eliminated entirely, because the 
H 3 f -02-S + (and H 3 f -02-C + ) cycles are fully quenched, and neutralization occurs entirely by 
grain recombination. 

These results strongly suggest that for realistic clouds, in which gas-grain neutralization 
is a likely important process, bistability will not occur. The precise value of k g in dark 
clouds may be variable, as the grain population can be modified by coagulation (Draine 
1985; Flower et al. 2005) and other processes. However, even with the removal of all grains 
with sizes below ~ 15 A, the rate coefficient k g ~ 10~ 15 cm 3 s~ 4 (e.g. McKee 1989), which is 
sufficiently large to suppress bistability. We thus confirm and explain the results of Shalabiea 
and Greenberg (1995), who found that bistability is eliminated when gas-grain processes, 
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including grain assisted recombination, are incorporated in the numerical computations. 

7. Summary 

The non-linearity of equations (l)-(4) is a mathematical prerequisite for bistability and 
the appearance of multiple solutions (Le Bourlot et al. 1993, 1995; Lee et al. 1998; Pineau 
des Forets & Roueff 2000; Charnley & Markwick 2003). However, the purpose of our paper 
has been to identify and analyze the chemical mechanisms underlying the bistability phe- 
nomenon. We have shown that a simple chemical instability, involving atomic ions (typically 
S + ), O2 molecules, and Hjj~ ions, allows multiple solutions to appear near the transition points 
from low-density high-ionization (HIP) to high-density low-ionization (LIP) conditions in 
the gas-phase chemistry. The instability drives the O2 density to either small (HIP) or large 
(LIP) values. Bistability occurs when a low O2 abundance in the HIP can be maintained 
up to maximal densities that are greater than the minimal densities above which a large 2 
density is maintained in the LIP. Two modes of neutralization are then simultaneously avail- 
able to the gas. The slow (HIP) mode is radiative recombination (usually S + + e — > S + u). 
The fast (LIP) mode is formation of molecular ions followed by dissociative recombination 
(usually S+ + 2 -> SO+ + O; SO+ + e S + O). 

We have derived simple analytic scaling relations for the range of gas densities and ion- 
ization rates for which bistability occurs. Our formulae show how the bistable range depends 
on the assumed gas-phase abundances, and several reaction rate-coefficients, including the 
dissociative recombination rate coefficient of H3 . Bistability is eliminated when the two 
modes of gas-phase neutralization are superseded by gas-grain processes such as grain as- 
sisted recombination of the atomic ions. This likely happens for the range of electron and 
O2 densities where bistability typically occurs in purely gas-phase models. We conclude that 
the bistability phenomenon probably does not occur in realistic dusty interstellar clouds, 
although it remains of theoretical interest. 
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n H2 /C_ 17 (cm 3 ) n H2 /C_ 17 (cm 3 ) 

Fig. 1. — The steady-state densities ni/C-17, and density fractions r^/n^, of electrons, S + , 
C + , and !Lt ions, C, and O atoms, and CO and O2 molecules, as functions of 71h 2 /C-i7- An 
HIP exists for densities less than 675 cm -3 , and a LIP exists for densities greater than 370 
cm -3 . The gas is bistable between these densities. 



-23- 




Fig. 2. — The H3 -C>2-S + cycle. Increasing the electron density, reduces H3 via dissociative 
recombination. This reduces O2 formed by the sequence (R6)-(R11). Destruction of S + 
by O2 decreases, so the S + density increases. This cycle can "run away" when S + is the 
dominant positive charge carrier (see text). 
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Fig. 3. — The fractional ionization, x e = n(e)/nn 2 , versus 71h 2 /C-17; f° r our standard model. 
The dominant destruction reactions for S + , , 2 , and C, are indicated for each of the 
three HIP, MIX, and LIP branches in the bistability region. E.g., "S + + e" indicates that 
S + is removed by recombination with electrons. The triangles indicate which "legs" of the 
H^-C^-S" 1- instability cycle operate in each of the three phases. 
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Fig. 4. — The fractional ionization, x e = n e /riu 2 , for five values of the H3 dissociative 
recombination rate coefficient ^4, ranging from 0.2 to 5 times our standard value of 1.7 x 10~ 7 
cm 3 s" 1 . For the lowest value of ^3,4 (dot-dashed curve) C + rather than S + is the dominant 
positive ion near the transition point in the HIP (see text). 
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Fig. 5. — The fractional ionization for grain-assisted recombination rate-coefficients for all 
atomic ions varying from kg = 1.0 x 1(T 17 to 1.0 x 10~ 14 cm s . For our standard model, 
k 9 = 0. 
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o 2 


4.2( 
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:-3) 
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:-3) 


0.55 


CN 


4.5( 
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:-6) 
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:-6) 
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2.2 


;-2) 
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:-2) 
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Note. — Densities (cm 3 ) of selected species, for the three HIP, MIX, and LIP solutions, 
for n H2 = 4.3 x 10 2 cm -3 , and ( = 1.0 x 1(T 17 s _1 . Column 5 lists HIP/LIP density ratios. 



